Nelkin scaling for the Burgers equation and the role of high-precision calculations 
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Nelkin scaling, the scaling of moments of velocity gradients in terms of the Reynolds number, is an 
alternative way of obtaining inertial-range information. It is shown numerically and theoretically for 
the Burgers equation that this procedure works already for Reynolds numbers of the order of 100 (or 
even lower when combined with a suitable extended self-similarity technique). At moderate Reynolds 
numbers, for the accurate determination of scaling exponents, it is crucial to use higher than double 
precision. Similar issues are likely to arise for three-dimensional Navier-Stokes simulations. 

PACS numbers: 47.27.-i,47.11.Kb,47.27. Jv 



Nelkinllfl, showed that the multifractal model of turbu- 
lence [2, 3], implies certain scaling relations for moments 
of velocity gradients (henceforth gradmoments) . Accord- 
ing to Nelkin, at high Reynolds numbers, when plotted as 
a function of the Reynolds number R, the pth moment of 
any component Vit of the velocity gradient should scale, 
to leading order, as 



((Vu) p ) 



(1) 



The exponents \v are expressible in terms of the multi- 
fractal structure function exponents £ p (cf. [l| or Sec. 
8.5.6). 

By using very highly resolved direct numerical simula- 
tion, it has been checked by Schumacher, Sreenivasan and 
Yakhot that not only is such scaling present (its first veri- 
fication), but that it is already seen at Reynolds numbers 
around 200, well below those where structure functions 
show any inertial-range scaling Q. This is perhaps not 
so suprising, given that inertial-range scaling is for in- 
termediate asymptotics with two large parameters, the 
Reynolds number and the ratio of the scale under con- 
sideration to the typical dissipation scale, whereas Nelkin 
scaling just requires a large Reynolds number. 

The one-dimensional Burgers equation 



dtu + ud x ii = vd x u; u(x, 0) = uq(x), 



(2) 



where u is the velocity and v the kinematic viscosity, 
can can throw light on why gradmoments display good 
scaling at such moderate Reynolds numbers. Further- 
more, it allows analytical determination of all the domi- 
nant and subdominant terms in the high-Reynolds num- 
ber expansion of gradmoments. We note that in a recent 
paper Q, the Burgers equation was used to illustrate 
why the extended self-similarity (ESS) technique Q gives 
improved scaling through the depletion of subdominant 
corrections. 

Heuristically, it is quite simple to show that for the 
Burgers equation we expect \p — P ~ !• Indeed, at high 
Reynolds numbers, the solutions of <j2j) display shocks 



broadened by viscosity over a distance 0{v) = O 
Within a shock, the pth power of the velocity gradient is 
O (R p ). Since shocks cover a fraction O (R 1 ) of the one- 
dimensional spatial domain, the stated scaling results. Of 
course such an argument tells us nothing about subdom- 
inant corrections and thus cannot be used to predict at 
what kind of Reynolds numbers this scaling emerges. 

We shall now address these issues more systematically, 
using simulations and theory. We shall also address a 
new question: scaling exponents are notoriously known 
with poor accuracy (cf., e.g., [3]); how accurately can 
we determine such exponents by working with Reynolds 
number at which there are significant subdominant 
corrections to scaling? Using recent results of van 
der Hoeven [1, 0], we shall show that this requires 
a subtle tradeoff between Reynolds numbers and pre- 
cision (number of decimal digits) used in the calculations. 

We begin with simulation-based results for the 
Reynolds number dependence of gradmoments when 
standard double-precision calculations suffice. We follow 
here the same strategy as in Ref. [6j : we solve the Burg- 
ers equation @ with the initial condition Uq(x) = sinx, 
using a pseudo-spectral method combined with fixed- 
time-step fourth-order Runge-Kutta time marching and 
a slaved scheme, known by the acronym ETDRK4 fioj ]. 
for handling the viscous dissipation. The gradmoments 
of integer order p, as a function of the Reynolds number 
R = 1/V, are defined as spatial averages over the period 
2tt: 



(3) 



Gradmoments are calculated for orders p from two to ten 
and Reynolds numbers R from twenty to one thousand. 
The number of collocation points N is taken between 8K 
and 256K where K stands for 2 10 = 1024; the time step St 
is between 10~ 5 and 10 -6 . We checked that the errors on 
gradmoments stemming from spatial and temporal trun- 
cation stay below the level needed for a double-precision 
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FIG. 1. (Colour online) Compensated pth-order moments (p 
from 2 to 10) of velocity gradient (M p ) versus both R (con- 
tinuous line with points in blue) and the ESS-type surrogate 
7? (dashed red line). 



calculation. The output is calculated at t = 2 when 
the solution of the Burgers equation has a well-developed 
shock. Since, as explained above, gradmoments are ex- 
pected to behave asymptotically as i? p_1 at large R, we 
display them in compensated manner, that is divide them 
by i? p_1 . Figure [1] shows the compensated gradmoments 
as a function of Reynolds number. Visual inspection 
shows that the expected flat behavior of the compen- 
sated gradmoments sets in around R = 40 for the lowest 
order p = 2 and around R = 300 for p — 10. In contrast, 
inertial-range scaling for structure functions, calculated 
from the same solution of the Burgers equation, appears 
clean only around Reynolds numbers of several thousands 
Q. This discrepancy, of nearly two orders of magnitude, 
can be made even larger by resorting to a procedure in- 
spired from ESS in which one resorts to a surrogate of 
the spatial separation, such as the third-order structure 
function and plots structure functions versus the surro- 
gate. In the case of gradmoments, we observe that the 
mean energy dissipation is given in terms of the mean 
square velocity gradient by e — vM 2 — (l/i?)M 2 . This 
has a finite positive limit Eoq as the Reynolds number 
tends to infinity. Hence, we can use R = Mijsao as a 
(suitably normalized) surrogate of the Reynolds number. 
This we call ESS-type plotting. The same Fig. [T] also 
shows this type of plotting. Now, the data look almost 
completely flat, except for the largest value of p around 
R = 20 where the data bend slightly upwards, as revealed 
by looking at the figure from the side [? ]. 

Of course, all this has to do with subdominant correc- 
tions to scaling and the way they arc affected by the ESS- 
type procedure. We now turn to theoretical interpreta- 
tions. For this we use the exact solution of the Burgers 
equation, obtained by employing the Hopf-Cole method 
[ill \v\ that transforms the Burgers equation into the 



heat equation. For the initial condition u(x, 0) = sin(x) 
this solution reads 

u(x,t) = -2vd x \n6(x,t) (4) 

p27T 

9(x,t)= c cos ^- x '^^G(x',t)dx'. (5) 
Jo 

Here, G(x',t) = Efc=-oc e ikx '~ uk2t is the Green's func- 
tion for the heat equation in the 27r- periodic case. We 
want to use this solution to determine the asymptotics of 
gradmoments for small v, i.e., large R. Using the method 
of steepest descent, in a way somewhat similar to what 
is found in Ref. fl3j, one can show that, for large R and 
any integer p > 2 

M P (R) = ApRP- 1 + B p R p - 2 + C p R p - 3 + ... (6) 

The coefficients are given by rather complicated and nu- 
merically ill-conditionned integrals. 

The expansion (J6j) and the numerical values of the co- 
efficients can actually be obtained by an alternative semi- 
numerical procedure, called asymptotic extrapolation, de- 
veloped by van der Hoeven [§] (see also Q for an elemen- 
tary presentation). Let us now say a few words about 
this technique, which will also be used below in connec- 
tion with high-precision spectral calculations. Suppose 
we have determined numerically with high precision the 
values of a function f(n) for integers n up to some high 
value N. We wish to obtain from this as many terms as 
possible in the high-n asymptotic expansion of /. Try- 
ing to fit the function by a guessed leading asymptotic 
form with some free parameters, will generally lead to 
very poor accuracy in such parameters. With some in- 
formation about the structure of the various terms in 
the expansion, a better method is to fit an expression 
containing one or several subdominant corrections (all 
with some unknown parameters) . Lacking such informa- 
tion, asymptotic extrapolation handles the problem by 
applying to the data a sequence of suitably chosen trans- 
formations that successively strip off the dominant and 
subdominant terms in the expansion for large n. At cer- 
tain stages of such transformations, the processed data 
allow simple extrapolations, most often by a constant. 
The transformations are meaningful as long as the succes- 
sively transformed data is free from conspicuous rounding 
noise and n has reached a simple asymptotic behavior 
(e.g. flat). From the extrapolation stages, it then be- 
comes possible (by undoing the transformations made) to 
obtain the asymptotic expansion of the data (including 
the values of the various parameters) up to some order 
which depends on the precision of the data and on the 
value of N. Here, we will denote the transformations by 
using the notation of Ref. Q. Thus, I stands for "in- 
verse", R for "ratio", SR for "second ratio" and D for 
"difference" . The sequence of transformations is chosen 
through various tests which provide some clue about the 
asymptotic class in which the data falls. 
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order (p) 


Xp 


Ap 


x ( p 1] 


Bp 


(2) 

Xp 


Cp 


2 


0.9999987 


+0.09032605 


- 0.002 


- 0.2290236 


- 1.002 


+0.2011 


3 


1.999998 


- 0.03245271 


1.00001 


+0.1736854 


0.005 


- 0.1325 


4 


2.999996 


+0.01249279 


2.00001 


- 0.090466 


1.0001 


+0.08417 


5 


3.999995 


- 0.00498725 


3.00001 


+0.045622 


1.99988 


- 0.08209 
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+ U.UD1UO 


7 


5.999993 


- 0.00084414 


5.000008 


+0.010955 


4.0002 


- 0.0398 


8 


6.999992 


+0.0003539 


5.999993 


- 0.00526 


5.002 


+0.024 


9 


7.999994 


- 0.0001495 


6.99991 


+0.0025 


6.009 


- 0.01 


10 


9.00001 


+0.000063 


7.9995 


- 0.0012 


7.03 


+0.03 



TABLE I. Dominant scaling exponents \p an d the first two subdominant exponents x\> an d Xp together with the corre- 
sponding coefficients A p , B p , and C p for the large- R behavior of gradmoments of order p, obtained by asymptotic extrapolation 
processing of a 400-digit precision determination of gradmoments from the Hopf-Cole solution. The theoretical values are 

Xp = P ~ 1. Xp = P ~ 2 . an <l Xp 2> = P - 3. 



To apply asymptotic extrapolation to the determina- 
tion of the coefficients in the high-Reynolds number ex- 
pansion (jGj) , we calculate the Hopf-Cole solution d4j)-d5j) 
and the gradmoments ([3]) using extreme precision float- 
ing point calculations [14| with 400 decimal digits. This 
precision guarantees that the only source of errors is lack 
of simple asymptoticity. The convolution structure of 
([5]) allows the use of fast Fourier transforms, also in very 
high precision [la ], for calculating 9, u and various space 
derivatives. The Reynolds number R is given all inte- 
ger values from 18 to i? ma x = 400. The processing of 
the gradmoments for p from 2 to 10 involves typically 
15 stages of transformations, the first eight of which are 
always R - 1, I, D, D, I, D, D, D [? ]. From the 
undoing of the transformations, using the "most asymp- 
totic data points" for determining constants, we obtain 
the following expansion: 

M P {R) = A P R X " + B p R x( p + C p R x( r +... (7) 

The results are shown in Table [U Only those digits of the 
coefficients that agree when processing the data succes- 
sively with i? max = 200 and R max — 400 are shown. It 
is seen that the scaling exponents for the dominant term 

X p and the first and second subdominant terms, Xp^ an d 
(2) 

Xp are very close to their theoretical values obtained 
from (j6)). The relative discrepancies are in the range 
10 -5 - 10 -6 for the dominant exponent and the accu- 
racy degrades for subdominant corrections, as expected. 

From the expansion (JT)) we can readily understand why 
Nelkin scaling appears at rather moderate Reynolds num- 
ber: the absolute value of relative correction stemming 
from the first subdominant term is R \B p / A p \. For ex- 
ample, it reaches the ten percent level which is easily 
picked up visually at R = 10\B p /A p \. Table [TTI shows 
the values of \B p /A p \ and we now understand why flat 
compensated gradmoments are seen in Fig. [1] beyond 
Reynolds numbers, varying with p, from a few tens to 
a few hundreds. To understand the ESS-type even bet- 
ter scaling, we expand the surrogate R in terms of R. 



From © with p — 2 and noticing that A2 — £00, we 
obtain 

R = R 1 + ^R° + 0(R- 1 ) . (8) 
A2 

Eliminating R between © and ([5]), we obtain 

Mp = ApRP- 1 + BpRP- 2 + ...-B p = Bp- ^~} )Ap B 2 . 

A2 

(9) 

Note that the expansion in terms of the surrogate R 
has the same structure as © and precisely the same 
dominant-term coefficient A p . However the coefficient B p 
of the first subdominant correction is significantly smaller 
than Bp (in absolute value) and may have a different sign. 
This explains for example why the graph for the compen- 
sated third-order gradmoment in terms of R bends down 
at the low end while it bends very slightly up in terms 
of R. As a consequence of the reduced subdominant cor- 
rections, the asymptotic behavior of gradmoments in the 
ESS-type representation emerges at Reynolds numbers 5 
to 20 times smaller than in the ordinary representation 
(see Table HI|. 

We should not be carried away and state that good 



order (p) 


r; = \B P /A P \ 


Rp — Bp 1 Ap | 


2 


2.5344 


0.0 


3 


5.3520 


0.2827 


4 


7.2414 


0.3622 


5 


9.1477 


0.9906 


6 


11.0613 


1.6116 


7 


12.9785 


2.2290 


8 


14.8980 


2.8440 


9 


16.8222 


3.4544 


10 


19.0604 


3.7507 



TABLE II. Estimates of Reynolds numbers beyond which sub- 
dominant corrections become small in the Reynolds number 
representation (middle column) and in the ESS-type repre- 
sentation (last column). 

scaling can emerge already at very moderate Reynolds 
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number provided we take the right quantity (here, grad- 
moments) and the right data processing technique (here, 
ESS). It all depends on what we call "good scaling". If 
we want to obtain scaling exponents with an error not ex- 
ceeding 10~ 2 or 10~ 3 , a flat looking compensated graph is 
definitely not enough since this is achieved as soon as the 
relative error is somewhere below 10 _1 . We now address 
the issue how asymptotic (how large in Reynolds num- 
ber) and how precise should a spectral calculation be in 
order to truly give accurate scaling exponents. Of course, 
the higher the Reynolds number, the lower the relative 
subdominant corrections will be. But, without enough 
precision, the simultaneous determination of dominant 
terms and subdominant corrections, say by asymptotic 
extrapolation, will be unable to handle more than very 
few such corrections and thus gives us substantial errors 
in the final results. In order to be closer to more realis- 
tic models such as the multi-dimensional Navier-Stokes 
equations, in investigating the trade-off between asymp- 
toticity and precision, we refrain from using the exact 
solution of the Burgers equation and resort to time in- 
tegration by (pseudo-)spectral technique. We use double 
and quadruple precision, both combined with asymptotic 
extrapolation, so as to obtain the most accurate possi- 
ble parameters. We calculate the scaling exponents %4 
and xe of the fourth and the sixth gradmoments, whose 
theoretical exact values are three and five, respectively. 
We determine how accurately we can predict these ex- 
ponents when applying asymptotic extrapolation (which 
for this purpose is substantially better than the aforemen- 
tioned ESS technique), using various maximum Reynolds 
numbers i? ma x- In double precision we were able to 
use three stages and in quadruple precision eight stages 
of the aforementioned transformations. The maximum 
wavenumber and the size of the time step are the same 
as reported at the beginning of the paper. We checked, by 
further halving of spatial and temporal resolutions, that 
they contribute negligible errors to the result. Figured] 
shows the relative errors for the two types of precision as 
function of i? ma x- It is striking that, when doubling the 
precision we can decrease the Reynolds number by about 
a factor of eleven (from 1 000 to 90) and still obtain a sub- 
stantial decrease (by a factor of 3 to 10) in the relative 
error. For accurate determination of scaling exponents, 
increasing the precision is here definitely more efficient 
than increasing the Reynolds number. It remains to be 
seen if this result carries over to a much broader class 
of equations, including multi-dimensional incompressible 
problems displaying random behavior. Already, we can 
state that the use of Nelkin scaling to analyze multifractal 
scaling in simulated 3D turbulent flow should definitely 
be encouraged, and preferably combined with high pre- 
cision caculations. 

We are indebted to Joris van der Hooven, T. Mat- 
sumuto, D. Mitra, O. Podvigina, and V. Zheligovsky for 
a number of useful discussions. S.C. thanks academic and 




FIG. 2. (Color online) Relative error of Nelkin exponents X4 
and X6 obtained by asymptotic extrapolation from pseudo- 
spectral calculations up to a maximum Reynolds number 
-Rmax- Upper set of curves: double precision calculations {xa- 
red filled circles, xs'- blue filled triangles); lower set of curves: 
quadruple precision (x<f- red inverted triangles, x&'- blue filled 
squares) . 



financial support rendered by NBIA (Copenhagen); and 
Danish Research Council for a FNU Grant No. 505100- 
50 - 30,168. The work was partially supported by ANR 
"OTARIE" BLAN07-2.183172. Some of the computa- 
tions used the Mesocentre de calcul of the Observatoire 
de la Cote d'Azur. 



* sagar@nbi.dk 
t uriel@oca.eu 

* Walter.Pauls@ds.mpg.de 

^ samriddhisankarray@gmail.com 

[1] M. Nelkin, Phys. Rev. A 42, 7226 (1990). 

[2] B. Mandelbrot, J. Fluid Mech. 62, 331 (1974). 

[3] G. Parisi and U. Frisch, in Turbulence and Predictabil- 
ity in Geophysical Fluid Dynamics and Climate Dynam- 
ics, edited by M. Ghil, R. Benzi and G. Parisi (North- 
Holland, Amsterdam, 1985), p. 84. 

[4] U. Frisch, Turbulence — The Legacy of A. N. Kol- 
mogorov, (Cambridge University Press, Cambridge, 
1995). 

[5] J. Schumacher, K. R. Sreenivasan and V. Yakhot, New 

J. Phys. 9, 89 (2007). 
[6] S. Chakraborty, U. Frisch and S. S. Ray, J. Fluid Mech. 

649, 275 (2010). 
[7] R. Benzi, S. Ciliberto, R. Tripiccione, C. Baudet, F. Mas- 

saioli and S. Succi, Phys. Rev. E 48, R29 (1993). 
[8] J. van der Hoeven, J. Symb. Comput. 44, 1000 (2009). 
[9] W. Pauls, and U. Frisch, J. Stat. Phys. 127, 1095 (2007). 



■5 



[10] S. M Cox and P. C. Matthews, J. Comp. Phys. 176, 430 
(2002). 

[11] E. Hopf, Comm. Pure Appl. Math. 3, 201 (1950). 
[12] J. D. Cole, Quart. Appl. Math. 9, 225 (1951). 
[13] D. O. Crighton and J. F. Scott, Phil. Trans. Roy. Soc. 
London 292, 101 (1979). 



[14] D. H. Bailey, "High-Precision Arithmetic in Scientific 
Computation" , Computing in Science and Engineer- 
ing, May- June, 2005, pg. 54-61; LBNL-57487. See also 
http://crd.lbl.gov/ ~dhbailey/ 

[15] http://www.kurims.kyoto-u.ac.jp/ ~ooura/fft.html. 



